27/11/20 PH
Data analysis for 450/440eV pump-probe, 1\mum focus (I think) - should be comparable to previous analysis posted by TMO team (Taran...?), running from v4 preprocessed data.
Using VMIproc class for VMI data analysis & interactive plotting, including uncertainties.
Implemented:
To do:
Background:
# # Memory profiler - OPTIONAL for testing
# # https://timothymonteath.com/articles/monitoring_memory_usage/
# %load_ext memory_profiler
# %memit
# # Quick hack for slow internet connection!
# %autosave 36000
# Standard imports
from pathlib import Path
import sys
import numpy as np
import xarray as xr
# Import class from file
sys.path.append(r'/cds/home/p/phockett/dev/tmo-dev/tmo')
import tmoDataBase as tb
# from vmi import VMI # VMI class
from inversion import VMIproc # VMI processing class
tb.setPlotDefaults() # Set plot defaults (optional)
# Set main data path
baseDir = Path('/reg/data/ana15/tmo/tmolw0618/scratch/preproc/v4')
# Setup class object and which runs to read.
data = VMIproc(fileBase = baseDir, runList = range(128,137+1), fileSchema='_v4')
# The class has a bunch of dictionaries holding info on the data
# data.runs['files']
# Read in the files with h5py
# There is very basic checking implemented currently, just to see if 'energies' is present.
data.readFiles()
# Add a filter to an existing filterset
# Set energy filter to match previous analysis
# TODO: implement event count filter here to check for saturation
# data.setFilter({'signal':{'energies':[0.05, 0.1]}})
data.setFilter(reset=True)
data.setFilter({'signal':{'evrs':[1,1,70]}}) # For v4 data, gas == evrs[:,70] true/false
data.setFilter({'bg':{'evrs':[0,0,70]}})
data.setFilter({'energies':[0.05, 0.1]})
data.filters
Note this defaults to electron coords ('xc','yc').
data.genVMIXmulti()
# Check an image with Xarray plotter (Matplotlib)
data.imgStack['signal'].sel(run=128).pipe(np.log10).plot.imshow()
# Check some metrics
print(data.data[129]['signal']['metrics'])
print(data.data[129]['bg']['metrics'])
# Subtraction, and set as new image set
data.imgStack['sub'] = data.imgStack['signal'] - data.imgStack['bg']
data.imgStack
# View full image stack
data.showImgSet()
This is set as native pixel value, or assumed to be the max point in the image.
# data.setCentre([553, 546]) # Set explictly
data.setCentre() # Default case - set as max ind.
data.checkCentre()
# All the coord parameters are set in `centreInds`
# Note here that `cList` gives the array indicies for the centre point, based on the pixel coords set.
data.centreInds
See the previous demo for info, here we'll just generate subtracted images.
# Subtraction, and set as new image set
data.imgStack['sub'] = data.imgStack['signal'] - data.imgStack['bg']
data.imgStack
Multiple methods to be supported... so far only Elio's cpBasex is implemented, and this needs a bit of setup.
# For cpbasex, set the basis path
basisPath = r'/reg/d/psdm/tmo/tmolw0618/scratch/results/tmo_ana/calc/G_r512_k128_l4.h5'
# Set also local version - pip installed version not working with loadG in testing (version mismatch issue, or something in import routine?)
pbasexPath = r'/reg/d/psdm/tmo/tmolw0618/results/modules/pbasex' #/pbasex'
# Import module + basis functions
data.setMethod(pbasexPath = pbasexPath, basisPath = basisPath)
Invert a set of images by setting filterSet (if not set, all images in self.imgStack will be processed). Optionally, set normalisation parameters for the images, the default is norm={'type':'max', 'scope':'global'}, which normalises to the global maximum (i.e. over all images).
Currently the options are:
type:
scope:
# Invert a dataset
data.inv(filterSet='signal')
# Full results are stored in a dictionay, self.proc[filterSet]
filterSet = 'signal'
data.proc[filterSet].keys()
# ... where 'xr' contains the results in an Xarray
data.proc[filterSet]['xr']
# For rapid plotting per run, use self.plotSpectra()
data.plotSpectra()
# This can also be set to return a HoloMap for more flexibility
hmap = data.plotSpectra(returnMap = True)
hmap.overlay('BLM').layout('run').cols(1)
Various other stages and results are also logged to the data strucure...
Images are set to 'filterSet-TYPE', where TYPE is:
quadrant.unfoldQuadrant for full symmetrized images.)These can be accessed and plotted with the usual functionality.
NOTE: there is currently no default radial masking set, so the default plot scaling will likely be off - adjust with the histogram.
# Quad filter
data.showImg(name='signal-fold')
data.showImg(name='signal-inv')
# Raw cpbasex output data
data.proc['signal']['pbasex'].keys()
A basis Poissonian bootstrapping routine is implemented for estimating uncertainties. This operates on the event data, and allows error propagation through the image generation & analysis, but is (consequently) a little slow. (For more on the method employed, see the wikipedeia page#Poisson_bootstrap).)
Note this routine is a little crude at the moment! It will run image generation for all filterSets, run image subtraction, then image inversion & stats for all filterSets, or just the listed element.
# Run bootstrapping & inversion routine, 5 iterations, and calculate stats for the subtracted data only.
data.bootstrapInv(N = 5, filterSet = 'sub')
# Stats are output to a dataset, with variables per iteration
data.stats
# Similarly to plotSpectra(), there is a plotUncertainties() function for this data
data.plotUncertainties()
# As before, an object can be returned for further plotting options
hmap = data.plotUncertainties(returnMap=True)
hmap.overlay('BLM') #.cols(1)
hmap.overlay('run').layout('BLM').cols(1)
# tmo-dev class version
data.__version__
import scooby
scooby.Report(additional=['holoviews', 'xarray'])